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ABSTRACT 

A van  der  Pol  relaxation  oscillator  that  is  subjected  to  external 
sinusoidal  forcing  can  exhibit  6table  and  unstable  periodic  and  almost 
periodic  responses.  For  some  forcing  amplitudes  it  even  happens  that 
two  stable  subharmonlcs  having  different  periods  may  coexist.  We  in- 
vestigate here  the  stable  responses  of  such  forced  oscillators.  By 
numerically  computing  the  rotation  number  of  stable  oscillations  for 
various  values  of  the  forcing  amplitude  and  oscillator  tuning,  we 
obtain  descriptions  of  regions  of  phase  locking,  successive  bifurcation 
of  stable  subharmonic  and  almost  periodic  oscillations,  and  overlap 
regions  where  two  distinct  stable  oscillations  can  coexist. 
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0.  Introduction. 

The  van  der  Pol  equation  is  typical  of  components  arising  in  a 
wide  variety  of  applications  ranging  from  electrical  circuitry  to  phy- 
siology. In  this  paper , we  study  it  in  the  fora 

2 2 2 

d u/dt  + k(u  - 1)  du/dt  + u B cos (^t  + cD 

where  k is  the  tuning  parameter,  and  B, p.  and  a are  the  (normalized) 
forcing  amplitude,  frequency  and  phase,  respectively. 

Periodic  solutions  of  this  equation,  and  their  stability  properties 
have  been  extensively  studied  in  the  past  fifty  years.  The  regular  case 
(k«l)  has  been  studied  with  averaging  methods  by  Bogoliuboff,  et.al.[l], 
and  by  analytic  and  topological  methods  by  Cartwright  and  Llttlewood  and 
others! 2, 3) . 

The  singular  case  (k»l)  is  more  difficult  to  analyze.  Early 
studies  of  the  free  problem  (B  - 0)  provided  descriptions  of  the  relax- 
ation oscillation  and  its  periodth].  The  period  of  the  free  oscillation 

was  shewn  to  be 

T ■ (3-2  log  2)  k + 7-01U  k‘l/3  - 1.325  + 0(log  k/  k) 

for  k »1.  Later  studies  gave  methods  for  approximating  the  relaxation 
oscillation  by  complicated  combinations  of  matched  asymptotic  expansions 
and  averaging  procedures  IU,5]. 

The  forced  singular  problem  (k  » 1,  B 4 0)  exhibits  remarkable 
solutions,  and  it  serves  as  the  main  point  of  the  investigations  by 
Cartwright  and  Llttlewood! 3 ,6-10) . A similar  problem,  but  with 
u2  - 1 replaced  by  slgnumtu2  - 1)  was  studied  by  Levinson  [11] 


by  patching  together  explicit  solutions.  These  investigations  showed  that 
there  are  parameter  values  among  0 < B < 2/3,  k»  1,  for  which  two 
stable  periodic  solutions  coexist.  These  are  distinct  subharmonics  having 
different  periods,  say  2n  + 1 times  the  forcing  period  27T/^c  . One 
striking  implication  of  this  was  pointed  out  by  Levinson:  For  fixed 

B,  k and  yw.  in  such  an  overlap  region,  there  are  two  stable  subharmonics 
having  periods  £q  - (2n  - l)2fC/y»A  and  £^  - (2n  + 1)2  , respectively. 

Given  any  sequence  of  zeros  and  ones,  let  marks  be  made  on  the  time  axis 
using  lengths  Jq  and  1^,  successively  as  prescribed  by  the  sequence.  Then 
there  is  a solute  ion  of  the  equation  having  zeros  near  and  only  near 
the  marks.  Thus,  there  is  a kind  of  randomness  exhibited  by  the  solution 

set.  These  solutions  will  be  discussed  later. 

\ 

Two  other  studies  of  the  forced  singular  case  that  are  related  to 
our  investigation  were  made  by  Hayashl[12]  and  in  [13]*  In  the  first, 
a numerical  calculation  was  performed  on  an  analog  device  for  fixed  k. 
Subharmonics  having  periods  2,3 ,k,5  and  6 times  the  forcing  period  were 
found  for  various  values  of  B and  f* , and  an  overlap  region  for  solutions 
having  periods  one  and  three  times  the  forcing  period  was  found.  In  the 
second,  a formal  construction  was  used  to  derive  necessary  conditions 
for  the  existence  of  certain  subharmonics.  It  was  shown  that  conditions 
for  various  subharmonics  could  be  satisfied  simultaneously  by  some  values 
of  B and  k. 

Our  methods  are  described  in  section  1 of  this  paper,  and  in  section  2, 
the  results  of  our  calculations  are  presented. 
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® AVAILABLE  copy 

1.  Methods. 

The  model  Is  studied  here  for  the  parameters 

satisfying 

a - 0,  /K  - 1,  0 < B < .8,  0 < 1/k  < .2. 

The  choice  a = 0 of  course  is  unimportant  since  this  can  always  be 
acomplished  by  a translation  of  the  time  variable.  Our  choice  of  * 1 
was  made  in  order  to  more  easily  relate  our  results  to  those  of  Cart- 
wright, Littlewood  and  Levinson.  The  calculations  carried  out  by  Hayashi 

suggest  that  similar  results  to  ours  would  be  obtained  for  any  fixed/r. 

v4i  Pel  e.yvit'a/ent 

Setting  £ = 1/k  and  integrating  the^  equation  once  leads  to  the^first 

order  system 

dy/dt  • (l/£)(y  - y3/3  - v) 

(1) 

dv/dt  « £y  - B cos  t 

It  is  this  system  that  we  study  here  for  0 < B < .8,  0 < £<  .2. 

Solutions  of  this  system  for  various  initial  data  (y(0),v(0))  define 
a transformation  of  the  (y,v)-plane  at  t * 0 to  the  (y,v)-plane  at  t * 2 ru, 

P : (y(0),v(0)) — »•  (y(2*)  ,v(2rr))  . 

This  is  the  Poincare  mapping,  and  invariant  sets  for  this  mapping  can 
describe  periodic  and  almost  periodic  solutions  of  the  equation.  For 
example,  a solution( subharmonic) having  period  3*27C  will  define  a three 
point  invariant  set  for  P.  If  this  subhartaonic  is  stable,  then  a 
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neighborhood  of  any  one  of  these  points  Mill  approach  the  point  under 
iterations  of  * P*P»P. 

A given  oscillatory  solution  can  be  described  further  by  its  rotation 
number.  This  is  defined  by 

f - lim  (l/2nW)  | arg(Pry  ) | , 

for  any  point g * (yo,vQ)  on  invariant  set  of  P corresponding 

to  the  solution.  The  notation  arg(yo,VQ)  denotes  the  positive  angle 
this  vector  makes  with  the  positive  y-axis.  For  example,  if  the  oscill- 
ation is  a subharmonic  having  period  3*2 <C,  then  its  rotation  number  can 
be  either  1/3  or  2/3.  In  general,  the  rotation  number  of  an  oscillatory 
solution  can  be  rational  or  irrational.  When^p  is  rational,  the  denomin- 
ator gives  the  period.  Also,  when p is  rational,  there  may  be  several 
distinct  solutions  having  that  rotation  number.  For  example,  if  £ - 1/3, 
there  is  a three  point  set  that  is  Invariant  under  the  Poincare  map, 
and  through  each  one  there  is  a solution  having  j>  = 1/3 • 

The  rotation  number  was  first  Introduced  by  Poincare,  and  it  has 
been  used  and  studied  in  many  contexts. since  then.  It  is  known  that 
the  rotation  number  of  a stable  oscillation  depends  continuously  on 
the  parameters  in  the  problem,  and  so  it  is  constant  over  a neighborhood 
of  parameter  values.  In  our  problem  there  are  both  stable  and  unstable 
oscillations  present  for  most  of  the  parameter  values.  Therefore,  the 
system  cannot  be  characterised  by  a single  rotation  number.  Even  when 
attention  is  restricted  to  stable  oscillations,  there  need  not  be  a 
unique  rotation  number  since  two  stable  oscillations  can  coexist. 
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Our  methods  consist  of  fixing  values  for  B and  £ , then  selecting  an 
initial  point  (y(0),v(0)).  Next,  we  solve  (1)  using  an  implementation  of 
Gear’ s code  due  to  Hindmarsh[ 15] • This  is  a variable  step,  variable  order 
numerical  Integrator.  When  this  solution  equilibrates (usually  to  six  or 
seven  place  accuracy  over  several  periods),  its  rotation  number  is  easily 
calculated.  Two  techniques  were  used  for  selecting  initial  data.  First, 
we  allow  £ and  B to  change  systematically  along  a straight  line  in  the 
(B, O-plane, with  each  new  calculation  using  the  last  value  from  the  previous 
one  as  an  initial  point.  In  this  way,  overlap  regions  were  determined 
by  passing  over  them  in  two  different  directions.  The  other  technique 
involved  fixing  g and  £,  and  then  calculating  the  rotation  number  for 
a sequence  of  initial  data  lying  near  the  relaxation  oscillation  of 
the  free  problem. 

2.  Numerical  evaluation  of^and  its  implications. 

Both  kinds  of  calculations  described  in  section  1 were  performed 
for  0 < B < .8  and  0 < £ < .2.  The  results  of  these  calculations  are 
summarized  in  Figure  1 by  a contour  mapping  of  as  a function  of  B and  £ . 

(FIGURE  1) 

The  free  problem  (B  * 0)  is  described  along  the  line  B ■ 0 in 
Figure  1.  Since  the  forcing  frequency  is  1 in  the  forced  problem,  the 
candidates  for  solutions  that  develop  into  subharmonics  are  those  free 
oscillations  whose  periods  are  Integer  multiples  of  2 7C.  These  occur 
at  p values  indicated  in  Table  1. 

(TABLE  1) 

These  are  also  indicated  on  Figure  1,  and  we  see  that  appropriate  sub- 
harmonics  are  observed  near  these  points  for  B near  zero.  The  values  in 
Tasble  1 cone  from  evaluation  of  the  first  few  terms  in  T given  in  section  0. 
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For  positive  values  of  B,  we  see  that  there  are  regions  over  which 
p is  constant.  The  profile  of  p for  B - .1,  .01  < £ < .1,  is  shown  in 
Figure  2. 

(FIGURE  2) 

This  indicates  that  f is  a monotonically  increasing  function  of  £ which 
is  constant  over  intervals.  Such  behavior  is  known  in  other  problems 
1 16] , and  it  illustrates  the  well  known  phenomenon  of  phase  locking,  where 
the  frequency  of  the  response  remains  constant  when  parameters  are  changed 
slightly. 

The  profile  of  £ for  B = .475  is  more  complicated  since  discontin- 
uities in  f occur.  Figure  3 indicates  that  is  continuous  for  £ > .04, 
but  that  it  is  double  valued  over  intervals  below  that. 

(FIGURE  3) 

The  overlap  Intervals  describe  parameter  values  where  there  is  not  a unique 
stable  response.  The  implications  of  overlaps  will  be  discussed  later. 

We  did  not  attempt  to  calculate  irrational  values  of  f.  However, 
they  presu*  ibly  occur,  and  they  correspond  to  almost  periodic  solutions 
of  the  equation  that  have  quite  small  parameter  sets  over  wich  they  are 
•table. 

An  indication  of  the  behavior  of  unstable  solutions  is  given  by 
the  behavior  of  stable  solutions  near  where  they  lose  stability.  In 
Figure  4,  five  examples  are  presented  which  show  this.  B is  fixed  at 
B ■ .475,  € is  decreased  through  the  region  where  J*  * 1/3.  A stable 

solution  having  period  3-2ffis  plotted  for  each  of  several  g values  by 
projecting  the  solution  back  onto  the  (y,v) -plane. 

(FIGURE  4) 

These  show  that  as  ( decreases,  the  solution  comes  successively  closer 
to  the  free  relaxation  oscillation,  at  the  some  time  layering  itself  over 
the  two  vertical  segments. 


3 . Summary 


The  main  motivation  for  this  investigation  is  the  study  of  stable 

£\»re 

oscillations  in  the  forced  van  der  Pol  equation. . this  is  a difficult 


problem  to  study  analytically,  we  have  calculated  the  rotation  numbers 


of  the  stable  responses  of  a van  der  Pol  relaxation  oscillator  to  sin' 


usoidal  forcing.  The  results  suggest  that  the  rotation  number  </> 
defines  a continuous,  but  piece-wise  constant  surface  except  in  overlap 


regions  where  it  is  double  valued,  having  what  resemble  folds 


ranges  where f is  single  valued  illustrate  the  phenomenon  of  phase  locking 
l.e.,  the  response  frequency  is  constant  over  a neighborhood  of  parameter 
values.  Various  phase  locking  regions  of  subharmonic  and  ultra-sub- 


harmonic responses  are  described  by  these  calculations 


Calculations  done  near  the  start  of  overlaps  show  that  the  surface 


steepens  as  the  overlap  is  approached.  While  the  overlap  regions  might 


be  described  by  simple  folds  in  the  surface  with  the  upper  and  lower 


branches  being  stable,  our  methods  give  no  indication  of  the  actual 


behavior  of  the  unstable  solutions  having  rotation  numbers  between  the 


two  stable  ones.  In  addition,  there  may  be  more  than  two  stable  oscilla 


tions  having  distinct  p values  coexisting,  although  we  were  able  to 
detect  only  two  stable  responses  in  any  overlap  region.  For  example 


stable  solutions  having  high  periods (>  100*2£  ) and  small  domains  of 


attraction  would  be  beyond  our  calculations.  Certainly,  if  the  form 


of  van  der  Pol's  equation  were  modified,  as  in  studies  of  multi-vibrators 


then  coexistence  of  more  than  two  stable  oscillations  would  be  expected 
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The  overlap  regions  are  quite  interesting.  Since  stable  oscillations 
having  different  frequencies  coexist  near  the  free  relaxation  oscillation, 
their  domains  of  attraction  must  be  intertwined  in  a very  complicated  way. 
In  particular,  the  system’s  behavior  is  very  sensitive  to  small  pertur- 
bations. Presumably  the  random  behavior  suggested  by  Levinson  is  to 
be  found  among  the  solution  set  having  initial  values  on  the  boundary  of 
the  domains  of  attraction. 

Finally,  our  calculations  are  consistent  with  the  results  established 
by  Cartwright  and  Littlewood.  Their  results  can  be  demonstrated  in 
Figure  1 by  fixing  £ at  a small  value,  say  .01,  and  then  observing  the 
f values  as  B changes.  They  show  that  for  sufficiently  small  g,  the 
interval  0 < B < .8  can  be  broken  up  into  three  disjoint  sets.  One  con- 
sists of  intervals  over  which  j>  is  single  valued,  another  consists  of 
Intervals  over  which p is  double  valued,  and  the  third  is  a small  set 
that  is  ignored.  The  intervals  of  single  and  double  values  of  f are 
clear  in  Figure  1,  and  the  remainder  consists  of  sets  where  takes  on 
irrational  values. 

From  another  point  of  view,  the  profile  of  the  f surface;|for  B 
fixed  and  g decreasing, can  be  interpreted  as  a bifurcation  diagram  where 
only  stable  branches  are  plotted.  For  example,  in  Figure  2 as  £ decreases 
from  0.1,  the  solution  having  period  3-2IT  loses  its  stablility  with  the 
appearence  of  other  oscillations.  Presumably  the  periodic  solution  is 
still  present  for  smaller  values  of  £,  but  it  is  no  longer  stable. 
Therefore,  if  Figure  2 were  viewed  as  being  a bifurcation  diagram,  then 
each  point  on  the  graph  should  be  extended  as  a straight  line  to  6=  0. 

If  the  unstable  branches  are  Included,  then  the  diagram  will  be  solid  above 
the  lowest  stable  value  since  periodic  solutions  presumably  persist  for 
decreasing  f. 
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Other  interesting  oscillators  can  be  constructed  using  information 
compliled  here.  For  example,  if  B is  allowed  to  vary  slowly  back  and 
foith  across  an  over  lap  region,  then  the  solution  would  have  one f value 
until  the  overlap  boundary  is  crossed,  when  it  would  equilibrate  rapidly 
to  a solution  having  the  other  Rvalue.  Then  as  B reverses,  the  opposite 
happens.  On  the  other  hand,  if  B and  £ vary  slowly  in  other  regions, 
the  solution  could  appear  to  be  chaotic. 

Some  indication  of  the  unstable  solutions  was  obtained  here  from 
calculated  solutions  for  parameter  values  near  where  they  lose  their 
stability.  This  indicates  that  the  solutions  are  layered  near  the  free 
relaxation  oscillation,  suggesting  an  onion  skin  configuration  reminiscent 
of  that  found  in  [14]  for  a two  dimensional  strange  attractor. 

We  have  not  attempted  in  this  paper  to  rigorously  derive  any  of 
the  results  suggested  by  our  calculations.  However,  the  calculations  do 
give  an  indication  of  the  nature  of  stable  responses  in  this  system, 
which  pose  difficult  analytic  problems. 


TABLE  1 


Protosubharmonics : 


to 


Free  oscillations (B  ■ 0)  having  period  21*  n 


n 

£ 

2 

.16 

3 

.095 

4 

.069 

5 

.054 

6 

.044 

7 

.037 

8 

.033 

9 

.029 

10 

.026 

11 

.024 

12 

.022 




Bibliography 


N.  Krylov#  N.  Bogoliubof f , Introduction  to  Nonlinear 
Mechanics,  Princeton  U.  Press,  1947  (Annals,  of  Math. 
Series,  #11). 

J.  J.  Stoker,  Nonlinear  Vibrations,  Interscience  Publ., 

New  York,  1950. 

J.  E.  Littlewood,  On  Nonlinear  Differential  Equations  of 
the  Second  Order  III,  Acta  Math.,  $7  (1957),  267-3 

W.  Wasow,  Asymptotic  Expansions  for  Ordinary  Differential 
Equations,  Interscience  Publ.,  New  York,  1965. 

* 

J.  D.  Cole,  Perturbation  Methods  in  Applied  Math.,  Ginn 
(Blaisdell) , Boston,  Mass.,  1968. 

M.  L.  Cartwright,  J.  E.  Littlewood,  Anns.  Math.,  5±  (1951), 
1-37. 

# Camb.  Phil.  Soc.  Proc.,  45  (1949),  495-501. 

4 

# Anns.  Math.#  48  (1947)#  472-494. 

% 

# J.  London  Math.  Soc.#  20  (1945),  180-189. 

M.  L.  Cartwright#  Forced  Oscillations  in  Nonlinear  Systems, 
Contrib.  to  the  Theory  of  Nonlinear  Oscillations,  Vol.  I, 
Princeton  U.  Press  (Ann.  of  Math.  Series,  #20). 

Levinson,  A Second  Order  Differential  Equation  with 
Singular  Solutions#  Anns.  Math.#  50  (1949),  127-152. 

C.  Hayashi,  Nonlinear  Oscillations  in  Physical  Systems# 
McGraw-Hill#  1964. 


4 


13.  J.  Grasman,  E.  J.  M.  Veling , G.  M.  Willems,  Relaxation 
Oscillations  Governed  by  a Van  der  Pol  Equation,  SIAM  J. 
Applied  Math,  3I_  (1976),  667-676. 

14.  M.  Henon,  A Two  Dimensional  Mapping  with  a Strange 
Attractor,  Springer  Lecture  Notes  in  Mathematics,  No.  565, 

15.  A.  C.  Hindmarsh,  Gear.  Ordinary  Differential  Equation 
Solver,  UCID-30001  (Rev.  3),  Lawrence  Livermore  Laboratory, 
December  1974,  Livermore,  CA  94550. 

•16.  J.  K.  Moser,  Stable  and  Random  Motions  in  Dynamical  Systems, 
Ann.  Math.  Studies,  Princeton  univ.  Press,  1973. 


v 

■ 

t 


I 


I 


5/17^ 

3/11 


#>  f>  ) sL 


0.75 


2.00  3.00 


EPS  = >14900  B = .47500 
ROTATION  NUMBER  =1./  3. 


EPS  = .12900  B = .47500 
ROTATION  HUMBER  =1./  3. 


ea*e  oo*2  00*1  oo*o  oo*i 


EPS  = .10900  B = .47S00 
ROTATION  NUMBER  =1./  3. 


00*6  00*2 


3.00 


2.00  3.00 


s , 


l 


\ 4 

EPS  = .04900  B = 
ROTATION  NUMBER  = 


+ 


•47500 
1./  3. 

* — - * 


Vc-i 


SECURITY  CLASSIFICATION  OF  THIS  PAGE  (Wtw%  Data  Entered) 


/^REPORT  DOCUMENTATION  PAGE 


#Jj  AFflfiR /-I  TR—T  7—"  j6fl  1 5 

j.  TITLE  (and  SuBnfi5^SS*^m^~m 


OAfiF  READ  INSTRUCTIONS 

BEFORE  COMPLETING  FORM 

2.  GOVT  ACCESSION  NO.  3-  RECIPIENT’S  CATALOG  NUMBER 


s-jyPEorsi 


[■REQUENCY  ^NTRAINMENT  OF  A FORCED  VAN  DER  / f 9)  Interim  /Udj-T  , , ) 

?OL  OSCILLATOR  „ * I U ^ * 1 

0 r~  • 1 — * » - — — — — — — j *.  performing  org.  import  NUMBER 


17.  AUTHORfA) 


^Oj Joseph  E. /Flaherty 
-"iF.  C. /Hoppensteadt 


— AFOSR-<<«i-2818  ~/73~ 


r».  performing  organization  name  and  address- 

Rensselaer  Polytechnic  Institute 
Department  of 'Mathematical  Sciences 
Troy,  New  York  12181 


H.  CONTROLLING  OFFICE  name  and  address 

Air  Force  Office  of  Scientific 
Research  (NM) , Bolling  AFB 
Washington,  DC  20332 


10.  PROGRAM  ELEMENT.  PROJECT.  TASK 
AREA  • WORK  UNIT  NUMBERS 

"V  61102F—  Z) 


QL. 


number  of  pages 

23 


MONITORING  AGENCY  NAME  A AOORESSf//  dllloront  from  Controlling  Olllco)  | IS.  SECURITY  CLASS,  (ol  (Jlta  raporlj 


UNCLASSIFIED 


»5a.  DECLASSIFICATION'' DOWNGRADING 
SCHEDULE 


6.  DISTRIBUTION  STATEMENT  (ol  Nil*  Ro port) 


Approved  for  public  release;  distribution  unlimited 


*7.  DISTRIBUTION  STATEMENT  (ol  Iho  otolroel  ontorod  In  Block  20.  II  dlllotont  from  Roport) 


•B.  SUPPLEMENTARY  NOTES 

Submitted  for  publication  in  the  Archive  for  Rational  Mechanics 
and  Analysis 


KEY  WORDS  (Continue  on  roeotoo  aldo  II  nocoaoory  and  Identity  by  block  number) 

forced  Van  der  Pol  oscillator,  relaxation  oscillations, 
frequency  entrainment,  phase  locking 


jected  to 
nstable 

rcing  amplitude 
different 
responses  of 
he  rotation 
the  forcing.—^  i 


DD  /jTn  1473 


UNCLASSIFIED 


SECURITY  CLASSIFICATION  or  This  PAGEfHTian  Data  Enltttd) 


amplitude  and  oscillator  tuning , we  obtain  descriptions; of 
regions  of  phase  locking,  successive  bifurcation  of  stable 
subharmonic  and  almost  periodic  osbilla'tions,  and  overlap 
regions  where  two  distinct  stable  oscillations  can  coexist. 
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